Code
library(zoib)
library(tidyverse)
library(GGally)
library(ggpubr)
library(kableExtra)
library(mice)
theme_set(theme_pubr(legend = "bottom"))
nnns_imputed<- readRDS("../Haojia_work/nnns_imputed.rds")Beta Regression
read imputed data
Code
dat=complete(nnns_imputed,1)
attach(dat);Percent_of_feeds_taken_by_mouth_at_discharge [1] 1.00 0.19 0.00 0.01 0.00 0.00 0.06 0.00 0.00 0.00 0.50 0.00 0.00 1.00 0.50
[16] 0.02 0.00 0.00 0.29 0.00 0.00 0.01 0.14 0.30 0.06 0.41 0.95 0.01 0.17 1.00
[31] 0.34 0.50 0.00 0.00 0.49 0.00 0.49 0.00 0.83 0.03 0.00 0.00 0.02 0.00 0.00
[46] 0.00 1.00 1.00 0.00 0.00 0.20 0.01 0.00 1.00 0.05 0.00 0.14 0.00 0.12 1.00
[61] 0.80 0.00 0.90 0.25 0.05 0.02 0.00 0.25 0.05 0.00 0.00 0.15 0.00 0.00 0.08
[76] 0.00 0.03 0.30 0.23 0.01 0.50 0.09 0.10 1.00 0.44 0.80 0.50 0.22 0.15 0.00
[91] 0.00 0.50 0.30 0.33 0.20 0.42 0.00 0.03 0.12 0.06 0.67 0.00 0.00 0.00 0.14
[106] 0.00 0.11 0.13 1.00 1.00 0.09 0.50 0.00 1.00 0.00Use ZOIB package
ZOIB Model
\[ f(y_i|\eta_i) = \begin{cases} p_i, & y_i =0 \\ (1-p_i)q_i, & y_i =1 \\ (1 − p_i)(1 − q_i)Beta(\alpha_{i_1}, \alpha_{i_2}) &y_i \in (0, 1) \end{cases} \]
\[ \begin{aligned} logit\left(\mu^{(0,1)} = \frac{\alpha_1}{\alpha_1+ \alpha_2}\right) = \mathbf x_1 \boldsymbol \beta_1 + I_1(\mathbf z_{1,i}\gamma_{1,i}) \\ log\left(v_i = \alpha_1+ \alpha_2\right) = \mathbf x_2 \boldsymbol \beta_2 + I_2(\mathbf z_{2,i}\gamma_{2,i}) \\ logit\left(p_i\right) = \mathbf x_3 \boldsymbol \beta_3 + I_1(\mathbf z_{3,i}\gamma_{3,i}) \\ logit\left(q_i\right) = \mathbf x_4 \boldsymbol \beta_4 + I_1(\mathbf z_{4,i}\gamma_{4 ,i}) \end{aligned} \]
Fit ZOIB model
Code
set.seed(11282023) #Set seed, bayesian model!
tictoc::tic()
pre_op_model =zoib(Percent_of_feeds_taken_by_mouth_at_discharge #Response
~Pre_Op_NNNS_attention_score+
Length_of_intubation_days+
Genetic_Syndrome_or_Chromosomal_Abnormality+
Cardiac_Anatomy+
Age_at_Surgery_days+
Female
#x1 design matrix
| 1 #x2 design matrix- estimating variance
|Pre_Op_NNNS_attention_score+
Length_of_intubation_days+
Genetic_Syndrome_or_Chromosomal_Abnormality+
Cardiac_Anatomy+
Age_at_Surgery_days+
Female #x3 design matrix
|Pre_Op_NNNS_attention_score+
Length_of_intubation_days+
Genetic_Syndrome_or_Chromosomal_Abnormality+
Cardiac_Anatomy+
Age_at_Surgery_days+
Female, #x4 design matrix
data = dat,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4, # number of MCMC chains
n.iter = 5000, #number of iterations before burn/thin
n.thin = 2, # thinning period
n.burn = 200, # burn-in period
seeds = c(11,29,20,23) #vector of seeds for chains to make reproducable
)
tictoc::toc()
saveRDS(pre_op_model, file="model1.RData")Interpreting output:
b: vector of estimates from Eqn 1; that is, g(mu) = xb*b +z*gamma
d: vector of estimates from Eqn 2; that is, log(eta) = xd*d+z*gamma
b0: vector of estimates from Eqn 3; that is, g(p0) = x0*b0 +z*gamma
b1: vector of estimates from Eqn ;4 that is, g(p1) = x1*b1+z*gamma
Code
model1 =readRDS(file="model1.RData")
model1$coeff |> summary()
Iterations = 1:2400
Thinning interval = 1
Number of chains = 4
Sample size per chain = 2400
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
b[1] -0.40816 0.74818 0.0076361 0.0075400
b[2] 0.08670 0.16883 0.0017231 0.0017843
b[3] -0.12256 0.06370 0.0006502 0.0007321
b[4] 0.01044 0.38861 0.0039662 0.0043701
b[5] 0.51766 0.39137 0.0039944 0.0040275
b[6] -0.22897 0.35027 0.0035749 0.0040599
b[7] -0.01677 0.03065 0.0003128 0.0003285
b[8] -0.44165 0.30822 0.0031458 0.0033937
b0[1] -1.88085 1.22359 0.0124882 0.0130445
b0[2] -0.14518 0.25755 0.0026286 0.0026731
b0[3] 0.27908 0.09623 0.0009821 0.0010380
b0[4] 0.31831 0.59714 0.0060946 0.0063665
b0[5] 1.54558 0.57721 0.0058911 0.0065491
b0[6] 1.07857 0.58879 0.0060093 0.0069353
b0[7] -0.04440 0.05236 0.0005344 0.0005509
b0[8] -0.70889 0.49425 0.0050444 0.0052464
b1[1] -5.26454 2.52180 0.0257380 0.0299346
b1[2] 0.60754 0.55635 0.0056782 0.0067350
b1[3] -0.38975 0.21021 0.0021454 0.0031057
b1[4] -72.91600 41.62390 0.4248222 7.9233266
b1[5] -0.09163 1.61224 0.0164549 0.0211416
b1[6] -0.49582 1.00420 0.0102491 0.0123597
b1[7] 0.24888 0.10494 0.0010711 0.0013713
b1[8] 1.47877 0.88376 0.0090198 0.0100330
d 1.01102 0.17214 0.0017569 0.0019064
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
b[1] -1.85699 -0.90675 -0.40844 0.087590 1.0715096
b[2] -0.23977 -0.02749 0.08513 0.200625 0.4193043
b[3] -0.25062 -0.16489 -0.12183 -0.079148 -0.0006274
b[4] -0.75772 -0.24505 0.01394 0.272360 0.7683315
b[5] -0.25352 0.26030 0.52092 0.780517 1.2806785
b[6] -0.91720 -0.46624 -0.22671 0.007560 0.4539170
b[7] -0.07966 -0.03713 -0.01625 0.003648 0.0429441
b[8] -1.05363 -0.64703 -0.44025 -0.230532 0.1575015
b0[1] -4.34565 -2.68888 -1.88189 -1.035950 0.4109481
b0[2] -0.65047 -0.31737 -0.14961 0.029522 0.3649150
b0[3] 0.10721 0.21125 0.27318 0.339820 0.4848661
b0[4] -0.84659 -0.07684 0.31247 0.717412 1.4897533
b0[5] 0.43129 1.15243 1.53832 1.928722 2.6996217
b0[6] -0.05431 0.68697 1.08013 1.465759 2.2500901
b0[7] -0.15423 -0.07724 -0.04288 -0.008099 0.0526682
b0[8] -1.68330 -1.04233 -0.70207 -0.367095 0.2396023
b1[1] -10.69062 -6.84610 -5.12750 -3.502697 -0.7335187
b1[2] -0.43126 0.23030 0.59253 0.959609 1.7767425
b1[3] -0.84252 -0.52501 -0.37801 -0.241454 -0.0087959
b1[4] -150.35996 -106.74572 -73.38067 -37.844099 -7.5182184
b1[5] -3.60955 -1.03421 0.03328 0.992953 2.7713079
b1[6] -2.49611 -1.15544 -0.49286 0.190776 1.4089553
b1[7] 0.06095 0.17607 0.24163 0.315919 0.4728664
b1[8] -0.22888 0.88836 1.45247 2.053894 3.2583483
d 0.66224 0.89529 1.01546 1.129489 1.3390093Check convergence
Traceplots
Code
traceplot(model1$coeff)
Autocorrelation plots
Code
autocorr.plot(model1$coeff)
